<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>clalsd.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="CLALSD.1"></a><a href="clalsd.f.html#CLALSD.1">CLALSD</a>( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
     $                   RANK, WORK, RWORK, IWORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      CHARACTER          UPLO
      INTEGER            INFO, LDB, N, NRHS, RANK, SMLSIZ
      REAL               RCOND
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      INTEGER            IWORK( * )
      REAL               D( * ), E( * ), RWORK( * )
      COMPLEX            B( LDB, * ), WORK( * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="CLALSD.22"></a><a href="clalsd.f.html#CLALSD.1">CLALSD</a> uses the singular value decomposition of A to solve the least
</span><span class="comment">*</span><span class="comment">  squares problem of finding X to minimize the Euclidean norm of each
</span><span class="comment">*</span><span class="comment">  column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
</span><span class="comment">*</span><span class="comment">  are N-by-NRHS. The solution X overwrites B.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The singular values of A smaller than RCOND times the largest
</span><span class="comment">*</span><span class="comment">  singular value are treated as zero in solving the least squares
</span><span class="comment">*</span><span class="comment">  problem; in this case a minimum norm solution is returned.
</span><span class="comment">*</span><span class="comment">  The actual singular values are returned in D in ascending order.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  This code makes very mild assumptions about floating point
</span><span class="comment">*</span><span class="comment">  arithmetic. It will work on machines with a guard digit in
</span><span class="comment">*</span><span class="comment">  add/subtract, or on those binary machines without guard digits
</span><span class="comment">*</span><span class="comment">  which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
</span><span class="comment">*</span><span class="comment">  It could conceivably fail on hexadecimal or decimal machines
</span><span class="comment">*</span><span class="comment">  without guard digits, but we know of none.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  UPLO   (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">         = 'U': D and E define an upper bidiagonal matrix.
</span><span class="comment">*</span><span class="comment">         = 'L': D and E define a  lower bidiagonal matrix.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  SMLSIZ (input) INTEGER
</span><span class="comment">*</span><span class="comment">         The maximum size of the subproblems at the bottom of the
</span><span class="comment">*</span><span class="comment">         computation tree.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N      (input) INTEGER
</span><span class="comment">*</span><span class="comment">         The dimension of the  bidiagonal matrix.  N &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  NRHS   (input) INTEGER
</span><span class="comment">*</span><span class="comment">         The number of columns of B. NRHS must be at least 1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  D      (input/output) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">         On entry D contains the main diagonal of the bidiagonal
</span><span class="comment">*</span><span class="comment">         matrix. On exit, if INFO = 0, D contains its singular values.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  E      (input/output) REAL array, dimension (N-1)
</span><span class="comment">*</span><span class="comment">         Contains the super-diagonal entries of the bidiagonal matrix.
</span><span class="comment">*</span><span class="comment">         On exit, E has been destroyed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  B      (input/output) COMPLEX array, dimension (LDB,NRHS)
</span><span class="comment">*</span><span class="comment">         On input, B contains the right hand sides of the least
</span><span class="comment">*</span><span class="comment">         squares problem. On output, B contains the solution X.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDB    (input) INTEGER
</span><span class="comment">*</span><span class="comment">         The leading dimension of B in the calling subprogram.
</span><span class="comment">*</span><span class="comment">         LDB must be at least max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RCOND  (input) REAL
</span><span class="comment">*</span><span class="comment">         The singular values of A less than or equal to RCOND times
</span><span class="comment">*</span><span class="comment">         the largest singular value are treated as zero in solving
</span><span class="comment">*</span><span class="comment">         the least squares problem. If RCOND is negative,
</span><span class="comment">*</span><span class="comment">         machine precision is used instead.
</span><span class="comment">*</span><span class="comment">         For example, if diag(S)*X=B were the least squares problem,
</span><span class="comment">*</span><span class="comment">         where diag(S) is a diagonal matrix of singular values, the
</span><span class="comment">*</span><span class="comment">         solution would be X(i) = B(i) / S(i) if S(i) is greater than
</span><span class="comment">*</span><span class="comment">         RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
</span><span class="comment">*</span><span class="comment">         RCOND*max(S).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RANK   (output) INTEGER
</span><span class="comment">*</span><span class="comment">         The number of singular values of A greater than RCOND times
</span><span class="comment">*</span><span class="comment">         the largest singular value.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK   (workspace) COMPLEX array, dimension (N * NRHS).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RWORK  (workspace) REAL array, dimension at least
</span><span class="comment">*</span><span class="comment">         (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS + (SMLSIZ+1)**2),
</span><span class="comment">*</span><span class="comment">         where
</span><span class="comment">*</span><span class="comment">         NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IWORK  (workspace) INTEGER array, dimension (3*N*NLVL + 11*N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO   (output) INTEGER
</span><span class="comment">*</span><span class="comment">         = 0:  successful exit.
</span><span class="comment">*</span><span class="comment">         &lt; 0:  if INFO = -i, the i-th argument had an illegal value.
</span><span class="comment">*</span><span class="comment">         &gt; 0:  The algorithm failed to compute an singular value while
</span><span class="comment">*</span><span class="comment">               working on the submatrix lying in rows and columns
</span><span class="comment">*</span><span class="comment">               INFO/(N+1) through MOD(INFO,N+1).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Based on contributions by
</span><span class="comment">*</span><span class="comment">     Ming Gu and Ren-Cang Li, Computer Science Division, University of
</span><span class="comment">*</span><span class="comment">       California at Berkeley, USA
</span><span class="comment">*</span><span class="comment">     Osni Marques, LBNL/NERSC, USA
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      REAL               ZERO, ONE, TWO
      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 )
      COMPLEX            CZERO
      PARAMETER          ( CZERO = ( 0.0E0, 0.0E0 ) )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      INTEGER            BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
     $                   GIVPTR, I, ICMPQ1, ICMPQ2, IRWB, IRWIB, IRWRB,
     $                   IRWU, IRWVT, IRWWRK, IWK, J, JCOL, JIMAG,
     $                   JREAL, JROW, K, NLVL, NM1, NRWORK, NSIZE, NSUB,
     $                   PERM, POLES, S, SIZEI, SMLSZP, SQRE, ST, ST1,
     $                   U, VT, Z
      REAL               CS, EPS, ORGNRM, R, RCND, SN, TOL
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      INTEGER            ISAMAX
      REAL               <a name="SLAMCH.130"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLANST.130"></a><a href="slanst.f.html#SLANST.1">SLANST</a>
      EXTERNAL           ISAMAX, <a name="SLAMCH.131"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>, <a name="SLANST.131"></a><a href="slanst.f.html#SLANST.1">SLANST</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           CCOPY, <a name="CLACPY.134"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>, <a name="CLALSA.134"></a><a href="clalsa.f.html#CLALSA.1">CLALSA</a>, <a name="CLASCL.134"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>, <a name="CLASET.134"></a><a href="claset.f.html#CLASET.1">CLASET</a>, <a name="CSROT.134"></a><a href="csrot.f.html#CSROT.1">CSROT</a>,
     $                   SGEMM, <a name="SLARTG.135"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>, <a name="SLASCL.135"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>, <a name="SLASDA.135"></a><a href="slasda.f.html#SLASDA.1">SLASDA</a>, <a name="SLASDQ.135"></a><a href="slasdq.f.html#SLASDQ.1">SLASDQ</a>, <a name="SLASET.135"></a><a href="slaset.f.html#SLASET.1">SLASET</a>,
     $                   <a name="SLASRT.136"></a><a href="slasrt.f.html#SLASRT.1">SLASRT</a>, <a name="XERBLA.136"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, AIMAG, CMPLX, INT, LOG, REAL, SIGN
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Test the input parameters.
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
<span class="comment">*</span><span class="comment">
</span>      IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( NRHS.LT.1 ) THEN
         INFO = -4
      ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN
         INFO = -8
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.155"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CLALSD.155"></a><a href="clalsd.f.html#CLALSD.1">CLALSD</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      EPS = <a name="SLAMCH.159"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'Epsilon'</span> )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set up the tolerance.
</span><span class="comment">*</span><span class="comment">
</span>      IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
         RCND = EPS
      ELSE
         RCND = RCOND
      END IF
<span class="comment">*</span><span class="comment">
</span>      RANK = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible.
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.0 ) THEN
         RETURN
      ELSE IF( N.EQ.1 ) THEN
         IF( D( 1 ).EQ.ZERO ) THEN
            CALL <a name="CLASET.177"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'A'</span>, 1, NRHS, CZERO, CZERO, B, LDB )
         ELSE
            RANK = 1
            CALL <a name="CLASCL.180"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>( <span class="string">'G'</span>, 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
            D( 1 ) = ABS( D( 1 ) )
         END IF
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Rotate the matrix if it is lower bidiagonal.
</span><span class="comment">*</span><span class="comment">
</span>      IF( UPLO.EQ.<span class="string">'L'</span> ) THEN
         DO 10 I = 1, N - 1
            CALL <a name="SLARTG.190"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( D( I ), E( I ), CS, SN, R )
            D( I ) = R
            E( I ) = SN*D( I+1 )
            D( I+1 ) = CS*D( I+1 )
            IF( NRHS.EQ.1 ) THEN
               CALL <a name="CSROT.195"></a><a href="csrot.f.html#CSROT.1">CSROT</a>( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
            ELSE
               RWORK( I*2-1 ) = CS
               RWORK( I*2 ) = SN
            END IF
   10    CONTINUE
         IF( NRHS.GT.1 ) THEN
            DO 30 I = 1, NRHS
               DO 20 J = 1, N - 1
                  CS = RWORK( J*2-1 )
                  SN = RWORK( J*2 )
                  CALL <a name="CSROT.206"></a><a href="csrot.f.html#CSROT.1">CSROT</a>( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
   20          CONTINUE
   30       CONTINUE
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Scale.
</span><span class="comment">*</span><span class="comment">
</span>      NM1 = N - 1
      ORGNRM = <a name="SLANST.215"></a><a href="slanst.f.html#SLANST.1">SLANST</a>( <span class="string">'M'</span>, N, D, E )
      IF( ORGNRM.EQ.ZERO ) THEN
         CALL <a name="CLASET.217"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'A'</span>, N, NRHS, CZERO, CZERO, B, LDB )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      CALL <a name="SLASCL.221"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
      CALL <a name="SLASCL.222"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If N is smaller than the minimum divide size SMLSIZ, then solve
</span><span class="comment">*</span><span class="comment">     the problem with another solver.
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.LE.SMLSIZ ) THEN
         IRWU = 1
         IRWVT = IRWU + N*N
         IRWWRK = IRWVT + N*N
         IRWRB = IRWWRK
         IRWIB = IRWRB + N*NRHS
         IRWB = IRWIB + N*NRHS
         CALL <a name="SLASET.234"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, RWORK( IRWU ), N )
         CALL <a name="SLASET.235"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, N, N, ZERO, ONE, RWORK( IRWVT ), N )
         CALL <a name="SLASDQ.236"></a><a href="slasdq.f.html#SLASDQ.1">SLASDQ</a>( <span class="string">'U'</span>, 0, N, N, N, 0, D, E, RWORK( IRWVT ), N,
     $                RWORK( IRWU ), N, RWORK( IRWWRK ), 1,
     $                RWORK( IRWWRK ), INFO )
         IF( INFO.NE.0 ) THEN
            RETURN
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        In the real version, B is passed to <a name="SLASDQ.243"></a><a href="slasdq.f.html#SLASDQ.1">SLASDQ</a> and multiplied
</span><span class="comment">*</span><span class="comment">        internally by Q'. Here B is complex and that product is
</span><span class="comment">*</span><span class="comment">        computed below in two steps (real and imaginary parts).
</span><span class="comment">*</span><span class="comment">
</span>         J = IRWB - 1
         DO 50 JCOL = 1, NRHS
            DO 40 JROW = 1, N
               J = J + 1
               RWORK( J ) = REAL( B( JROW, JCOL ) )
   40       CONTINUE
   50    CONTINUE
         CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, N, NRHS, N, ONE, RWORK( IRWU ), N,
     $               RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
         J = IRWB - 1
         DO 70 JCOL = 1, NRHS
            DO 60 JROW = 1, N
               J = J + 1
               RWORK( J ) = AIMAG( B( JROW, JCOL ) )
   60       CONTINUE
   70    CONTINUE
         CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, N, NRHS, N, ONE, RWORK( IRWU ), N,
     $               RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
         JREAL = IRWRB - 1
         JIMAG = IRWIB - 1
         DO 90 JCOL = 1, NRHS
            DO 80 JROW = 1, N
               JREAL = JREAL + 1
               JIMAG = JIMAG + 1
               B( JROW, JCOL ) = CMPLX( RWORK( JREAL ), RWORK( JIMAG ) )
   80       CONTINUE
   90    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         TOL = RCND*ABS( D( ISAMAX( N, D, 1 ) ) )
         DO 100 I = 1, N
            IF( D( I ).LE.TOL ) THEN
               CALL <a name="CLASET.278"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'A'</span>, 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
            ELSE
               CALL <a name="CLASCL.280"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>( <span class="string">'G'</span>, 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
     $                      LDB, INFO )
               RANK = RANK + 1
            END IF
  100    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Since B is complex, the following call to SGEMM is performed
</span><span class="comment">*</span><span class="comment">        in two steps (real and imaginary parts). That is for V * B
</span><span class="comment">*</span><span class="comment">        (in the real version of the code V' is stored in WORK).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
</span><span class="comment">*</span><span class="comment">    $               WORK( NWORK ), N )
</span><span class="comment">*</span><span class="comment">
</span>         J = IRWB - 1
         DO 120 JCOL = 1, NRHS
            DO 110 JROW = 1, N
               J = J + 1
               RWORK( J ) = REAL( B( JROW, JCOL ) )
  110       CONTINUE
  120    CONTINUE
         CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, N, NRHS, N, ONE, RWORK( IRWVT ), N,
     $               RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
         J = IRWB - 1
         DO 140 JCOL = 1, NRHS
            DO 130 JROW = 1, N
               J = J + 1
               RWORK( J ) = AIMAG( B( JROW, JCOL ) )
  130       CONTINUE
  140    CONTINUE
         CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, N, NRHS, N, ONE, RWORK( IRWVT ), N,
     $               RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
         JREAL = IRWRB - 1
         JIMAG = IRWIB - 1
         DO 160 JCOL = 1, NRHS
            DO 150 JROW = 1, N
               JREAL = JREAL + 1
               JIMAG = JIMAG + 1
               B( JROW, JCOL ) = CMPLX( RWORK( JREAL ), RWORK( JIMAG ) )
  150       CONTINUE
  160    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Unscale.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLASCL.323"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
         CALL <a name="SLASRT.324"></a><a href="slasrt.f.html#SLASRT.1">SLASRT</a>( <span class="string">'D'</span>, N, D, INFO )
         CALL <a name="CLASCL.325"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>( <span class="string">'G'</span>, 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
<span class="comment">*</span><span class="comment">
</span>         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Book-keeping and setting up some constants.
</span><span class="comment">*</span><span class="comment">
</span>      NLVL = INT( LOG( REAL( N ) / REAL( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
<span class="comment">*</span><span class="comment">
</span>      SMLSZP = SMLSIZ + 1
<span class="comment">*</span><span class="comment">
</span>      U = 1
      VT = 1 + SMLSIZ*N
      DIFL = VT + SMLSZP*N
      DIFR = DIFL + NLVL*N
      Z = DIFR + NLVL*N*2
      C = Z + NLVL*N
      S = C + N
      POLES = S + N
      GIVNUM = POLES + 2*NLVL*N
      NRWORK = GIVNUM + 2*NLVL*N
      BX = 1
<span class="comment">*</span><span class="comment">
</span>      IRWRB = NRWORK
      IRWIB = IRWRB + SMLSIZ*NRHS
      IRWB = IRWIB + SMLSIZ*NRHS
<span class="comment">*</span><span class="comment">
</span>      SIZEI = 1 + N
      K = SIZEI + N
      GIVPTR = K + N
      PERM = GIVPTR + N
      GIVCOL = PERM + NLVL*N
      IWK = GIVCOL + NLVL*N*2
<span class="comment">*</span><span class="comment">
</span>      ST = 1
      SQRE = 0
      ICMPQ1 = 1
      ICMPQ2 = 0
      NSUB = 0
<span class="comment">*</span><span class="comment">
</span>      DO 170 I = 1, N
         IF( ABS( D( I ) ).LT.EPS ) THEN
            D( I ) = SIGN( EPS, D( I ) )
         END IF
  170 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      DO 240 I = 1, NM1
         IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
            NSUB = NSUB + 1
            IWORK( NSUB ) = ST
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Subproblem found. First determine its size and then
</span><span class="comment">*</span><span class="comment">           apply divide and conquer on it.
</span><span class="comment">*</span><span class="comment">
</span>            IF( I.LT.NM1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              A subproblem with E(I) small for I &lt; NM1.
</span><span class="comment">*</span><span class="comment">
</span>               NSIZE = I - ST + 1
               IWORK( SIZEI+NSUB-1 ) = NSIZE
            ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              A subproblem with E(NM1) not too small but I = NM1.
</span><span class="comment">*</span><span class="comment">
</span>               NSIZE = N - ST + 1
               IWORK( SIZEI+NSUB-1 ) = NSIZE
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              A subproblem with E(NM1) small. This implies an
</span><span class="comment">*</span><span class="comment">              1-by-1 subproblem at D(N), which is not solved
</span><span class="comment">*</span><span class="comment">              explicitly.
</span><span class="comment">*</span><span class="comment">
</span>               NSIZE = I - ST + 1
               IWORK( SIZEI+NSUB-1 ) = NSIZE
               NSUB = NSUB + 1
               IWORK( NSUB ) = N
               IWORK( SIZEI+NSUB-1 ) = 1
               CALL CCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N )
            END IF
            ST1 = ST - 1
            IF( NSIZE.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              This is a 1-by-1 subproblem and is not solved
</span><span class="comment">*</span><span class="comment">              explicitly.
</span><span class="comment">*</span><span class="comment">
</span>               CALL CCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N )
            ELSE IF( NSIZE.LE.SMLSIZ ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              This is a small subproblem and is solved by <a name="SLASDQ.413"></a><a href="slasdq.f.html#SLASDQ.1">SLASDQ</a>.
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="SLASET.415"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, NSIZE, NSIZE, ZERO, ONE,
     $                      RWORK( VT+ST1 ), N )
               CALL <a name="SLASET.417"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'A'</span>, NSIZE, NSIZE, ZERO, ONE,
     $                      RWORK( U+ST1 ), N )
               CALL <a name="SLASDQ.419"></a><a href="slasdq.f.html#SLASDQ.1">SLASDQ</a>( <span class="string">'U'</span>, 0, NSIZE, NSIZE, NSIZE, 0, D( ST ),
     $                      E( ST ), RWORK( VT+ST1 ), N, RWORK( U+ST1 ),
     $                      N, RWORK( NRWORK ), 1, RWORK( NRWORK ),
     $                      INFO )
               IF( INFO.NE.0 ) THEN
                  RETURN
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              In the real version, B is passed to <a name="SLASDQ.427"></a><a href="slasdq.f.html#SLASDQ.1">SLASDQ</a> and multiplied
</span><span class="comment">*</span><span class="comment">              internally by Q'. Here B is complex and that product is
</span><span class="comment">*</span><span class="comment">              computed below in two steps (real and imaginary parts).
</span><span class="comment">*</span><span class="comment">
</span>               J = IRWB - 1
               DO 190 JCOL = 1, NRHS
                  DO 180 JROW = ST, ST + NSIZE - 1
                     J = J + 1
                     RWORK( J ) = REAL( B( JROW, JCOL ) )
  180             CONTINUE
  190          CONTINUE
               CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NSIZE, NRHS, NSIZE, ONE,
     $                     RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
     $                     ZERO, RWORK( IRWRB ), NSIZE )
               J = IRWB - 1
               DO 210 JCOL = 1, NRHS
                  DO 200 JROW = ST, ST + NSIZE - 1
                     J = J + 1
                     RWORK( J ) = AIMAG( B( JROW, JCOL ) )
  200             CONTINUE
  210          CONTINUE
               CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NSIZE, NRHS, NSIZE, ONE,
     $                     RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
     $                     ZERO, RWORK( IRWIB ), NSIZE )
               JREAL = IRWRB - 1
               JIMAG = IRWIB - 1
               DO 230 JCOL = 1, NRHS
                  DO 220 JROW = ST, ST + NSIZE - 1
                     JREAL = JREAL + 1
                     JIMAG = JIMAG + 1
                     B( JROW, JCOL ) = CMPLX( RWORK( JREAL ),
     $                                 RWORK( JIMAG ) )
  220             CONTINUE
  230          CONTINUE
<span class="comment">*</span><span class="comment">
</span>               CALL <a name="CLACPY.462"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'A'</span>, NSIZE, NRHS, B( ST, 1 ), LDB,
     $                      WORK( BX+ST1 ), N )
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              A large problem. Solve it using divide and conquer.
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="SLASDA.468"></a><a href="slasda.f.html#SLASDA.1">SLASDA</a>( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ),
     $                      E( ST ), RWORK( U+ST1 ), N, RWORK( VT+ST1 ),
     $                      IWORK( K+ST1 ), RWORK( DIFL+ST1 ),
     $                      RWORK( DIFR+ST1 ), RWORK( Z+ST1 ),
     $                      RWORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ),
     $                      IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ),
     $                      RWORK( GIVNUM+ST1 ), RWORK( C+ST1 ),
     $                      RWORK( S+ST1 ), RWORK( NRWORK ),
     $                      IWORK( IWK ), INFO )
               IF( INFO.NE.0 ) THEN
                  RETURN
               END IF
               BXST = BX + ST1
               CALL <a name="CLALSA.481"></a><a href="clalsa.f.html#CLALSA.1">CLALSA</a>( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ),
     $                      LDB, WORK( BXST ), N, RWORK( U+ST1 ), N,
     $                      RWORK( VT+ST1 ), IWORK( K+ST1 ),
     $                      RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
     $                      RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
     $                      IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
     $                      IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
     $                      RWORK( C+ST1 ), RWORK( S+ST1 ),
     $                      RWORK( NRWORK ), IWORK( IWK ), INFO )
               IF( INFO.NE.0 ) THEN
                  RETURN
               END IF
            END IF
            ST = I + 1
         END IF
  240 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Apply the singular values and treat the tiny ones as zero.
</span><span class="comment">*</span><span class="comment">
</span>      TOL = RCND*ABS( D( ISAMAX( N, D, 1 ) ) )
<span class="comment">*</span><span class="comment">
</span>      DO 250 I = 1, N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Some of the elements in D can be negative because 1-by-1
</span><span class="comment">*</span><span class="comment">        subproblems were not solved explicitly.
</span><span class="comment">*</span><span class="comment">
</span>         IF( ABS( D( I ) ).LE.TOL ) THEN
            CALL <a name="CLASET.508"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'A'</span>, 1, NRHS, CZERO, CZERO, WORK( BX+I-1 ), N )
         ELSE
            RANK = RANK + 1
            CALL <a name="CLASCL.511"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>( <span class="string">'G'</span>, 0, 0, D( I ), ONE, 1, NRHS,
     $                   WORK( BX+I-1 ), N, INFO )
         END IF
         D( I ) = ABS( D( I ) )
  250 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Now apply back the right singular vectors.
</span><span class="comment">*</span><span class="comment">
</span>      ICMPQ2 = 1
      DO 320 I = 1, NSUB
         ST = IWORK( I )
         ST1 = ST - 1
         NSIZE = IWORK( SIZEI+I-1 )
         BXST = BX + ST1
         IF( NSIZE.EQ.1 ) THEN
            CALL CCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB )
         ELSE IF( NSIZE.LE.SMLSIZ ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Since B and BX are complex, the following call to SGEMM
</span><span class="comment">*</span><span class="comment">           is performed in two steps (real and imaginary parts).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
</span><span class="comment">*</span><span class="comment">    $                  RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO,
</span><span class="comment">*</span><span class="comment">    $                  B( ST, 1 ), LDB )
</span><span class="comment">*</span><span class="comment">
</span>            J = BXST - N - 1
            JREAL = IRWB - 1
            DO 270 JCOL = 1, NRHS
               J = J + N
               DO 260 JROW = 1, NSIZE
                  JREAL = JREAL + 1
                  RWORK( JREAL ) = REAL( WORK( J+JROW ) )
  260          CONTINUE
  270       CONTINUE
            CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NSIZE, NRHS, NSIZE, ONE,
     $                  RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
     $                  RWORK( IRWRB ), NSIZE )
            J = BXST - N - 1
            JIMAG = IRWB - 1
            DO 290 JCOL = 1, NRHS
               J = J + N
               DO 280 JROW = 1, NSIZE
                  JIMAG = JIMAG + 1
                  RWORK( JIMAG ) = AIMAG( WORK( J+JROW ) )
  280          CONTINUE
  290       CONTINUE
            CALL SGEMM( <span class="string">'T'</span>, <span class="string">'N'</span>, NSIZE, NRHS, NSIZE, ONE,
     $                  RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
     $                  RWORK( IRWIB ), NSIZE )
            JREAL = IRWRB - 1
            JIMAG = IRWIB - 1
            DO 310 JCOL = 1, NRHS
               DO 300 JROW = ST, ST + NSIZE - 1
                  JREAL = JREAL + 1
                  JIMAG = JIMAG + 1
                  B( JROW, JCOL ) = CMPLX( RWORK( JREAL ),
     $                              RWORK( JIMAG ) )
  300          CONTINUE
  310       CONTINUE
         ELSE
            CALL <a name="CLALSA.571"></a><a href="clalsa.f.html#CLALSA.1">CLALSA</a>( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N,
     $                   B( ST, 1 ), LDB, RWORK( U+ST1 ), N,
     $                   RWORK( VT+ST1 ), IWORK( K+ST1 ),
     $                   RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
     $                   RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
     $                   IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
     $                   IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
     $                   RWORK( C+ST1 ), RWORK( S+ST1 ),
     $                   RWORK( NRWORK ), IWORK( IWK ), INFO )
            IF( INFO.NE.0 ) THEN
               RETURN
            END IF
         END IF
  320 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Unscale and sort the singular values.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="SLASCL.588"></a><a href="slascl.f.html#SLASCL.1">SLASCL</a>( <span class="string">'G'</span>, 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
      CALL <a name="SLASRT.589"></a><a href="slasrt.f.html#SLASRT.1">SLASRT</a>( <span class="string">'D'</span>, N, D, INFO )
      CALL <a name="CLASCL.590"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>( <span class="string">'G'</span>, 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="CLALSD.594"></a><a href="clalsd.f.html#CLALSD.1">CLALSD</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
